Measuring Estuarine Total Exchange Flow From Discrete Observations

Abstract The exchange between estuaries and the coastal ocean is a key dynamical driver impacting nutrient and phytoplankton concentrations and regulating estuarine residence time, hypoxia, and acidification. Estuarine exchange flows can be particularly challenging to monitor because many systems have strong vertical and lateral velocity shear and sharp gradients in water properties that vary over space and time, requiring high‐resolution measurements in order to accurately constrain the flux. The total exchange flow (TEF) method provides detailed information about the salinity structure of the exchange, but requires observations (or model resolution) that resolve the time and spatial co‐variability of salinity and currents. The goal of this analysis is to provide recommendations for measuring TEF with the most efficient spatial sampling resolution. Results from three realistic hydrodynamic models were investigated. These model domains included three estuary types: a bay (San Diego Bay), a salt‐wedge (Columbia River), and a fjord (Salish Sea). Model fields were sampled using three different mooring strategies, varying the number of mooring locations (lateral resolution) and sample depths (vertical resolution) with each method. The exchange volume transport was more sensitive than salinity to the sampling resolution. Most (>90%) of the exchange flow magnitude was captured by three to four moorings evenly distributed across the estuarine channel with a minimum threshold of 1–5 sample depths, which varied depending on the vertical stratification. These results can improve our ability to observe and monitor the exchange and transport of water masses efficiently with limited resources.

TEF has been widely applied in estuarine research and exchange flows more generally. The TEF framework has been used to determine freshwater fluxes from a small groundwater-driven estuary (Ganju et al., 2012), to estimate estuarine residence times (Lemagie & Lerczak, 2014;Sutherland et al., 2011), to examine the relationship between exchange flow and mixing (Wang et al., 2017), and study seasonal variability (Conroy et al., 2020;Giddings & MacCready, 2017), among others. Most of the aforementioned examples are modeling studies, with the exception of Ganju et al. (2012). There are also analyses of salinity flux from observations that do not use the TEF framework (e.g., Lerczak et al., 2006;MacDonald & Horner-Devine, 2008), but calculations of salt flux from observations are limited due to the large data requirement to resolve the temporal and spatial co-variability in salinity and currents as well as a lack of knowledge regarding flux errors when undersampling occurs.
The goal here is to provide recommendations for applying TEF to in situ observations, specifically to understand the most efficient spatial sampling resolution and the percent of the exchange flow captured under various strategies. This paper examines TEF calculated from sub-sampling realistic numerical models, representative of moorings in a channel, compared to TEF calculated from the full model resolution in order to compare how quickly the two estimates converge as the number of moorings increase. Three estuaries were studied in order to span much of the parameter space of estuarine characteristics. The objectives of this study were (a) to test how TEF converged for different sampling resolutions; (b) to examine how this varied between estuaries and sampling strategies; (c) to attempt to outline best practices for how many moorings and instruments would be required to quantify TEF from observations; and (d) to understand flux errors (magnitude and potential bias) when a cross-section is under-sampled. The three realistic estuary models, details of the TEF calculation, and the sampling methods are described in the methods Section 2. The current and salinity patterns characteristic of each estuary and individual cross-section are included in Section 3. The rest of the results are organized into sections based on the various sampling approaches: evenly distributed moorings (Section 4), strategically distributed moorings (Section 5), and a case study designed to approximate a simple observational approach (Section 6).

Realistic Hydrodynamic Models
Realistic hydrodynamic numerical models of three estuaries and their adjacent coastal regions were used ( Figure 1). These span different estuary types and geometries and include a small bay (San Diego Bay), salt-wedge 10.1029/2022JC018960 3 of 21 (Columbia River), and large fjord (Strait of Juan de Fuca in the Salish Sea; e.g., Geyer & MacCready, 2014). Results were extracted hourly from two across-channel sections in each model over a full year of simulation time in order to resolve the tides and capture seasonal variation. Extracted data is available on-line (Lemagie et al., 2022). Further details about each model are outlined in the following paragraphs. The Salish Sea model including the Strait of Juan de Fuca employed the Regional Ocean Modeling System (ROMS, Shchepetkin & McWilliams, 2005). Simulations from 2004 to 2007 were developed by the University of Washington Coastal Modeling Group MacCready et al., 2009;Sutherland et al., 2011). The model was forced with realistic river flow, tides, wind stress, surface heat flux, and open boundary conditions (e.g., Giddings et al., 2014) with initial and open boundary values for tracers, subtidal velocity, and subtidal surface height from the Navy Coastal Ocean Model (NCOM) Kara et al., 2006). The domain spans the inland waters of the Salish Sea (including Puget Sound, the Strait of Georgia, and the Strait of Juan de Fuca) and coastal ocean from 43N to 50N and 200 km offshore with a horizontal resolution of 1.5 km at the coast to 4.5 km far offshore. There were 40 sigma layers with enhanced vertical resolution near the surface and bottom. This analysis focuses on data extracted from 2005 at two cross-sections spanning the Strait of Juan de Fuca (SJF): near the ocean entrance (SJF EH1 ), and 16 km upstream (SJF J2C ; Figure 1). These sections correspond with previous analyses and validation of the model results . At SJF EH1 and SJF J2C the channel is 22.1 and 21.5 km wide, respectively, and 73 and 60 m deep (Table 1). Model skill was high (≥0.92) relative to observed currents, tidal sea surface elevation, salinity and temperature, although overall slightly too salty and cold (by ∼1.5 psu and ∼0.5°C) within the Salish Sea . Most pertinent to this study, the exchange flow through the Strait of Juan de Fuca compared well with observations and was insensitive to model resolution (Giddings & MacCready, 2017).   Barron et al., 2006). The domain extended from 39N to 50N and ∼300 km in the offshore direction with horizontal resolution from tens of meters in the estuary and river to 3 km in the ocean (Karna & Baptista, 2016;Karna et al., 2015). The vertical grid consisted of 37 sigma levels between sea level and 100 m datum and an equipotential z-grid below 100 m. Data were extracted from two cross-sections: the river mouth (CR M ) and 14 km upstream at site Saturn-03 (CR Sat03 ), to match previous studies where model validation was performed (Figure 1; e.g., Karna et al., 2015). At CR M and CR Sat03 the channel is 4.2 and 5.7 km wide and 18 and 16 m deep, respectively (Table 1). The model demonstrated high skill compared to long term observations, particularly outside of high discharge and neap tide conditions which are estimated to occur only 16% of the time (Karna & Baptista, 2016).
The San Diego Bay (SDB) and adjacent coastal dynamics were simulated using the Coupled Ocea n-Atmosphere-Wave-Sediment-Transport (COAWST) model system (Kumar et al., 2012;Warner et al., 2010) to represent the surfzone and shelf circulation (Wu et al., 2020(Wu et al., , 2021. This model grid sits within three one-way nested parent models using ROMS (Shchepetkin & McWilliams, 2005) and is coupled with the Simulating Waves Nearshore (SWAN) model to include surface gravity waves (Booij et al., 1999). Boundary and initial conditions for the outermost domain were from the California State Estimate (CASE) solution (Marshall et al., 1997) with tides from the Advanced CIRculation tidal database (Westernik et al., 1993) and surface forcing from the North American Mesoscale Forecast (NAM) and the Coupled Ocean-Atmosphere Mesoscale System (COAMPS). The largest grid extends between 29N and 36N and over 500 km offshore with 2 km horizontal resolution, which is downscaled to the finest grid with horizontal resolution from 8 m near the coast to 110 m at the western boundary, and has 10 stretched vertical sigma levels (Wu et al., 2020(Wu et al., , 2021. This study focuses on 2016 results at the estuary mouth (SDB M ) and 4 km upstream (SDB C ; Figure 1). At SDB M and SDB C , the channel is 0.5 and 1.2 km wide and up to 19 and 14 m deep, respectively (Table 1). This model has not been rigorously validated against observations of San Diego Bay, but exhibits circulation similar to prior observations .

Estuarine State Estimates
The vertical stratification index was used to characterize the water column at each cross section by where the overbar denotes a vertical mean. Vertical mean density was computed   Simpson et al. (1981). ϕ gives an estimate of the potential energy of the water column relative to the mixed state such that in a vertically-well mixed water column ϕ = 0. The influence of salinity and temperature on the density structure were computed by substituting = ( ( ), ) and = ( ( ) ) , respectively. While available potential energy is another useful framework for understanding estuarine systems (e.g., MacCready & Giddings, 2016), ϕ is useful in this context because as ϕ approaches 0 S in and S out converge. ϕ was computed for each lateral column separately before calculating an area-weighted cross-sectional mean value.
The Kelvin number K e and Ekman number E k provide an estimate of the degree of horizontal and vertical variability in the currents (Valle-Levinson, 2008) and may help predict how many moorings and vertical sample depths are needed to accurately constrain the exchange. The Kelvin number estimates the importance of Earth's rotation on the flow. Wide basins (K e > 2) are more likely to have strong horizontal shear (Garvine, 1995). The Ekman number estimates the importance of vertical mixing (Kasai et al., 2000;Winant, 2004). Large Ekman number (E k > 1) basins are likely to have strong horizontal shear regardless of their width (Valle-Levinson, 2008). The Kelvin and Ekman numbers were calculated from the full resolution of model fields and then time-averaged. = ( ′ ) − 1 2 for estuary width W, depth H, reduced gravity g′, and coriolis parameter f.
where A z is the flow's eddy viscosity. The eddy viscosity was not available from the Columbia River SELFE model, and E k could not be explicitly calculated.

TEF Calculations
Subtidal exchange flow is calculated using isohaline coordinates following the TEF dividing salinity method (Lorenz et al., 2019). Following the TEF framework (MacCready, 2011), the net transport of a tracer c through a cross-sectional area A(S > S′) determined by salinity S′ is defined as: where u is the velocity normal to the cross-section (positive values are into the estuary), t is time, and 〈〉 denotes a subtidal filter (here the Godin low pass filter, Thomson & Emery, 2014). A profile of the tracer exchange can also be determined by differentiation: The transport profile was separated into distinct inflow and outflow layers (l) by finding the extrema in the Q c profiles, ignoring extrema below a certain noise threshold Q thresh (Lorenz et al., 2019). Q thresh was defined here as a fixed percentage of the maximum transport magnitude, Q thresh = Q percent * max(|Q c (S)|), where Q percent = 0.01.
The salinity values associated with the Q c extrema-along with the salinity endpoints S min , S max -made up the dividing salinities S div and the transport in each layer was Inflow was positive and layers were defined by the sign of the net transport: Subscripts a and b are used to enumerate inflow and outflow layers respectively, following Lorenz et al. (2019).
The net exchanges were defined by: This summation transforms the results from l = 1: L individual layers into two layers, which does not greatly impact the result if the flow is predominantly two-layered. The mean inflow and outflow salinities can be calculated by: where the tracer, c, is the salinity, S and no superscript implies the volume flux only, ( ) = ⟨ ∫ ( ′ ) ⟩ . The above exchange flows and corresponding salinities are referred to as the TEF bulk values (e.g., Lorenz et al., 2019).
Currents and salinities were extracted hourly from each cross-section. Since the SELFE model uses an unstructured grid, CR output were interpolated onto horizontally fixed straight cross-sections at CR Sat03 and CR M , roughly matching the spatial resolution of the model grid. ROMS variables were extracted at the grid resolution. Velocities were rotated onto along-and across-channel coordinates, defined by the angle of the cross-section. The principle axis of the area-averaged currents over the year were closely aligned with each cross-section (e.g., Table 2). In order to avoid tidal aliasing, the start and end times were estimated by the timing of the spring tidal sea level maximum closest to each calendar end point.

TEF Sampling Strategies
Four methods of sampling the cross-sectional fields and calculating TEF were compared: (a) using horizontal and vertical resolution from the IxJ model grid to calculate TEF IJ ; (b) using an M × N array of evenly spaced samples to calculate TEF MN from an increasing integer number of "moorings" M evenly distributed across the channel width with N sample depths each, which were evenly distributed across the time-averaged channel depth at each location x = x m (e.g., Figure 2a); (c) using μ moorings with strategic placement of each mooring μ = 1, …, I determined by maximizing the correlation between TEF μJ and TEF IJ ; and (d) a case study TEF case designed to imitate observations with a single bottom-mounted Acoustic Doppler Current Profiler (ADCP) and salinity measurements near the surface and bottom of the water column at M evenly distributed mooring locations (e.g., Figure 2b). In cases 2-4, the width of the channel represented by each mooring is defined by the channel boundaries and the mid-point between adjacent moorings ( Figure 2). In cases 2 and 4, mooring locations shift as M varies so that all moorings are evenly spaced. Each of these methods is described in more detail in the following paragraphs.
One goal of this analysis was to identify the minimum number of moorings and samples required for various sampling approaches to converge to TEF IJ . An appropriate definition of convergence between these methods could depend on the specific application or research question. One comparison that demonstrated utility was to identify the threshold for which the magnitude of each TEF parameter consistently remained within ≤10% of the TEF IJ value. However, since the observed salinity range was small relative to the magnitude of salinity values ( Note. Along-channel flow was defined as positive flowing into the estuary normal to the cross-section and the principle axis is reported here as degrees counter-clockwise from section normal (with the normal vector directed into the estuary).

Table 2 Oceanographic Characteristics, Including the Time and Area-Averaged Mean Salinity, and Along-Channel Current Magnitude, the Principle Axis of the Currents, and the Kelvin and Ekman Numbers Calculated at Each Cross Section From Unfiltered Time Series
fraction. S in,out were converted to equivalent freshwater fractions using the tidal maximum salinity over time at each cross-section, S max following FWFin,out = ( max − in,out ) ( max) −1 .

Full TEF, TEF IJ
TEF IJ used the full vertical and horizontal resolution from the model grid. TEF IJ represented the expected value which other estimates of TEF were hypothesized to converge toward at high sampling resolutions. Velocity and salinity fields were sampled at points (x i , z ij (t)) and the area represented by each sample was computed by ΔA i,j (t) = Δx i Δz ij (t). Subscripts i and j indicate the indices on the model grid in the across-channel and vertical direction, respectively (e.g., Figure 2).

Evenly Distributed Subsamples, TEF MN
TEF MN used an MxN array of samples evenly distributed across the channel and throughout the water column (e.g., Figure 2a). This method was chosen to test the convergence of TEF MN toward TEF IJ as the number of moorings (M ≤ I) or sample depths (N ≤ J) were increased. This method is simple enough to be consistently applied in every case. However, when there is sharp bathymetric variability within width Δx m , it is not obvious how to estimate area ΔA m,n (t). For example, on a steep slope a gridded area could either over-estimate or under-estimate the flux. To address this, two approaches were compared: the first method assumed that u, S were constant with depth over distance Δx m and the second assumed the profiles of u, S had a consistent shape over distance Δx m and were thus constant across σ-levels as in Lerczak et al. (2006). The difference in the results between approaches was negligible. The results reported herein assume that u and S were constant along σ-levels to estimate ΔA MN (Equation 3) as illustrated on Figure 2a.

Strategically Located Subsamples, TEF μJ
TEF μJ was calculated by incrementally adding moorings in order based on identifying the mooring which contributed the largest improvement in the correlation coefficient between TEF μJ and TEF IJ , similar to the approach of using maximum explained variance used by Wei et al. (2020). Lateral mooring placement was sampled at μ ≤ I grid locations while J sampling depths were included at each mooring location. Since TEF is a derived flux quantity, it was necessary to estimate the cross-sectional area represented by each sample of u and S. This computation used linear interpolation assuming u and S were constant across σ-levels, analogous to TEF MN . The lateral edges of the regions represented by each mooring μ were defined by the channel edges and the mid-point between each mooring pair in the across-channel direction. Importantly, the moorings are not necessarily evenly spaced in this approach. The sample area ΔA μJ was calculated as the total model grid area at a given σ-level between the lateral edges bounding each mooring μ.

Case Study, TEF case
A case study was designed to imitate a sampling plan where S and u observations are not co-located and are constrained by common instrument and deployment logistics (Figure 2b). M moorings were evenly distributed laterally over the cross-section. At each mooring u was sampled at the full model resolution, to approximate having a bottom-mounted ADCP and salinity was sampled 1 m off of the bottom and 1 m below sea level, to approximate having a bottom-mounted sensor as well as one mounted from a surface float. From this sampling distribution two variations of vertical salinity interpolation were compared. In case A, S was linearly interpolated to the velocity sample depths (i.e., the model grid). In case B, a two-layer system was assumed having well-mixed surface and bottom layers each with constant S. The depth of the boundary between the well-mixed layers was approximated by the mean depth of the 0-crossing between inflow and outflow in the deepest part of the channel (50, 8, and 5 m at sections SJF EH1 , CR M , and SDB M , respectively). With observations this interpolation could be calculated during the analysis stage using observed currents, therefore this estimate of the mixed layer depth does not rely on a priori knowledge. In reality, vertical patterns of currents and salinities may be decoupled, or may vary over time, however this simplified approach is applied for the case study to be most relevant to an observational study with limited to no a priori knowledge about the system. While not always dynamically appropriate, it is a reasonable simplified approach for many estuarine systems (e.g., Aristizábal & Chant, 2015;Lerczak et al., 2006). Tests with extrapolated currents in the top and bottom 10% of the water column, and 2-5 m from the bottom-to simulate ADCP limitations-had negligible impact on the results. Additional case studies mimicking other sampling strategies were not included since TEF MN and TEF case already span the parameter space for most sampling approaches and some approaches, such as shipboard transects, are typically limited in duration.

Discrete Calculations
For the discrete calculation of Equation 3 applied to each TEF method the spatial coordinates were first converted to isohaline coordinates. The salinity range was defined by the minimum and maximum salinity sampled at each cross-section over the year. This range was divided into N bins = 500 evenly spaced salinity bins. The currents u, salinity S, and area A were interpolated onto the spatial grid defined for each method (e.g., Figure 2) and then mapped into these discrete salinity bins at each time prior to the calculation of TEF (Equation 3).

Exchange Flow
The estuaries and individual cross-sections chosen for this study differ in the degree of stratification and shear, the range of seasonal and tidal variability, as well as in the channel width and bathymetric complexity. These features contribute to differences in TEF. Before presenting a comparison of the TEF calculated by different methods and resolutions, the characteristics of the salinity, along-channel currents, and TEF at the full model resolution is discussed.

Strait of Juan de Fuca
Annual mean currents and salinity in the Strait of Juan de Fuca generally exhibit a classical pattern of estuarine circulation with outflow and relatively fresher water near the surface as well as inflow and saltier water at depth (Figure 3a). Occasional intrusions of the Columbia River plume during prolonged downwelling-favorable winds (Giddings & MacCready, 2017;Hickey et al., 2009) were apparent in the mean currents as an upstream flow near the surface at SJF EH1 ( Figure 3a); there was also a small, intermittent intrusion of fresher coastal water (<30 psu), associated with increased horizontal shear (Giddings & MacCready, 2017). This fresh surface inflow was not evident in the mean currents at SJF J2C (Figure 3b), which is further from the mouth. Neither the salinity (Figure 4a) nor stratification (Figure 4b) had strong subtidal variability, although temperature had a greater contribution to the vertical stratification in the latter months of the year than in the early spring. Annual mean TEF was predominantly two-layered, with outflow of fresher water and inflow of relatively saltier water (Figures 5d and 5e) with an occasional third inflow layer at SJF EH1 (Figures 5a and 5d).

Columbia River
At CR M near the Columbia River mouth there was also a mean outflow of fresher water near the surface and saltier near-bottom inflow (Figure 3c). The bathymetry is more complicated at CR Sat03 where a shallow sill bisects the channel. At CR Sat03 the mean inflow is weak and the mean currents are mostly out of the estuary (Figure 3d). At CR M there was a seasonal cycle in the salinity with a half-amplitude of 5 psu, similar to the mean tidal range (the salinity time series was low-pass filtered using a Godin filter; the scale of tidal variability is indicated in red). The stratification at CR M had large variability (similar to the mean) at seasonal, spring-neap, and higher tidal frequencies (Figure 4d). Note that stratification is calculated from salinity only. The salinity range of the inflow S in was greater at CR Sat03 than at CR M , with more evenly distributed transport across the salinity range, while the outflow across CR Sat03 was predominantly at salinities <10 psu. TEF was two-layered with outflow of fresher water and inflow of relatively saltier water (Figures 5f and 5g). At both sections the outflow Q out had greater seasonal variability than the inflow Q in (Figure 5b).

San Diego Bay
In San Diego Bay, the annual mean exchange was out of the estuary near the surface, but the surface waters are saltier than at the bottom (Figures 3e and 3f). In this shallow system with relatively little rainfall, vertical stratification at the mouth is thermally controlled much of the year  and varies at seasonal and spring-neap tidal frequencies (e.g., Figure 4f). Also, the spatial standard deviation of salinity across the section SDB M was small relative to the fluctuation of the mean. Variability in this system is driven by surface heating and evaporation as well as tidal advection and river discharge (Largier, 1995;Largier et al., 1997). TEF was only examined in isohaline coordinates for this study although an analysis in joint temperature-salinity space is also possible and may be valuable for a dynamical study (Lorenz et al., 2020). The seasonal variability in the exchange at SDB M near the mouth was reflected in the results. In winter the inflow salinity S in was <0.3 psu saltier than S out (Figure 5h). When the stratification was thermally dominated, the isohaline coordinate TEF reversed and the inflow was fresher than the outflow, again by a small margin. This pattern was similar further upstream at SDB C , but at SDB C the magnitude of the exchange, was weaker and to identify three distinct layers of the exchange flow as shown in Figure 5i, the threshold used to calculate the dividing salinities was adjusted to Q percent = 0.1, instead of 0.01. For consistency in the following sections, the same value of Q percent = 0.01 was used for all of the simulations. The time series of TEF parameters Q in and S in (Q out and S out ) were summed (and averaged, following Equation 8) across all inflow (outflow) layers and were not sensitive to Q thresh . Dashed lines mark the dividing salinities between inflow and outflow layers and colored dots indicate S in (red and orange to distinguish inflow layers) and S out (blue and cyan to distinguish outflow layers) values. Note that the exchange threshold used to calculate the transport in each layer following Equation 3 is Q thresh = 0.01, except as shown for cross-section SDB C (i), where for the same threshold there were 18 dividing salinities. For SDB C shown here Q thresh is 0.1.

Horizontal and Vertical Shear
The degree of horizontal shear relative to vertical shear is likely important for estimating how many moorings and sample depths would be needed to capture the exchange flow. At SJF EH1 and SJF J2C the channel is wide enough that the Coriolis force can be important (K e ∼ 2, Table 2) as can be seen in the tilted isopycnals at SJF J2C (Figure 3b). A detailed analysis of the mechanisms driving TEF at SJF EH1 and SJF J2C found that the flow was in geostrophic balance, but also that temporal changes were driven by the baroclinic pressure gradient and advection (Giddings & MacCready, 2017). In the Columbia River, the low Ke < 2 suggests that horizontal shear is likely to be small relative to vertical shear (Valle-Levinson, 2008). However, at CR Sat03 the steep bathymetry that divides the flow between two channels contributes to horizontal shear. While the role of friction is expected to be stronger in the shallow San Diego Bay, E k approaches the moderate frictional regime where both horizontal and vertical shear can be found (Valle-Levinson, 2008). Channel curvature (Figure 1) may also contribute to the observed across-channel shear and salinity gradients (Figures 3e and 3f). Overall, all of the cross-sections presented here exhibit both horizontal and vertical shear.

Evenly Distributed Moorings, TEF MN
TEF MN parameters Q in,out converged toward TEF IJ parameters Q in,out as the number of moorings M increased as long as there was a minimum number of sample depths N at each mooring ( Figure 6). The minimum threshold of sample depths for convergence between TEF MN and TEF IJ was N ≤ 4 in most cases (Table 3). Sampling the model fields at fewer depths resulted in smaller Q in,out magnitudes. For a small number of sample depths, for example, N < 4, TEF MN diverged from TEF IJ and approached 0 as the number of moorings M increased. At sections SJF EH1 and SDB M Q in,out depended on whether the channel center was sampled; the magnitudes were overestimated if only a single centered mooring was sampled and were underestimated if only two moorings were sampled (i.e., not sampling the channel center; Figure 6). Similarly, there was a lower correlation between Q in,out from TEF IJ and TEF MN for M = 2 than for M = 1 (Figure 7).
The deviations of S in,out from TEF IJ values were relatively small ( Figure 6)-within 10% of the freshwater fraction even with only a single mooring, in some cases (Table 3). Particularly at SJF EH1 and SJF J2C , where there is strong stratification, there was a greater range in S in,out as the number of depths N varied than as the number of moorings M varied.
If too few depths were sampled the annual mean magnitude of Q in,out and S in,out calculated from TEF MN did not converge to TEF IJ , but the number of sample depths required for convergence did not change as M increased (e.g., Figure 6a). In most cases the TEF MN parameters S in,out converged toward TEF IJ parameters at values of M and N that were similar to or smaller than the number of moorings and sample depths over which Q in,out converged. Thus, Q in,out were the limiting parameters with the TEF MN method.

Strategically Distributed Moorings, TEF μJ
In order to assess the sensitivity of TEF on sampling method, a strategic approach to lateral mooring placement rather than a geometric distribution was tested. The strategic approach, TEF μJ , incrementally added moorings that contributed the maximum correlation improvement between the sampled parameters, TEF μJ and TEF IJ . Due to the iterative nature of this method, only variability in the lateral mooring placement was assessed, and the water column was sampled at the full vertical model resolution, N = J. While the maximum correlation for a given number of moorings was slightly less with fewer vertical samples (N < J), the patterns were similar to those presented here for N = J, particularly for N ≥ 4.
Using the time series correlation to strategically select lateral mooring locations, the correlation between the sub-sampled TEF μJ and TEF IJ parameters converged to 1 for fewer moorings than TEF MN parameters (Figure 7). However, this distinction is minimal for the cross-sections in the Columbia River and San Diego Bay (e.g., Figures 7e and 7l) where the correlation is high (>0.8) even when only one or two moorings were used to sample the exchange flow. The high correlation between TEF IJ and that from a single optimal mooring (i.e., TEF μ = 1,J ) in the CR and SDB is likely due to the temporal variability in the exchange flow having a wide spatial signal, that is, similar temporal variability over the full cross-section. This is opposed to SJF where temporal variations in exchange have a strong spatial signature, such as those associated with intermittent downwelling-favorable winds (e.g., Figure 3a and Giddings & MacCready, 2017).
The mooring order varied based on which TEF parameter they were designed to capture (i.e., Q in , Q out , S in or S out ; Figure 8). However, the placement of the first mooring tended to be in deeper parts of the channel. In general the strategic mooring placement spanned the section. If each cross-section was geometrically divided into thirds, the first triad of moorings (i.e., the three tallest bars of each color on Figure 8) was roughly distributed across those three sections, resulting in a sampling distribution that was similar between TEF μ = 3,J and TEF M = 3,N . Despite this similarity, and despite the correlation converging more quickly for strategic moorings (Figure 7), the minimum number of moorings μ before Q in,out converged to within 10% of TEF IJ values was greater for the strategic sampling approach (Table 4) than when using the evenly spaced sampling approach (Table 3). This implies that the areas of high TEF variance do not fully correspond with the maximum TEF magnitude.
This investigation of TEF μJ provided insight into the sensitivity of the results to the specific sample locations by comparing the variation in results using the mooring order from each TEF μJ parameter (Figure 9). The results presented in Figures 7 and 8 and throughout the text primarily focus on the outcomes of each parameter of TEF μJ with strategic moorings selected based on that same parameter. This means that the single mooring (μ = 1) with the highest correlation to TEF IJ for the parameter Q in is not necessarily the same mooring location with the highest correlation to TEF IJ for the parameters Q out , S in , or S out . Similarly, the second mooring (μ = 2) is not necessarily the same across parameters Q in , Q out , S in , or S out and so on. For example, there was some variation in the pattern of Q in for increasing μ, between the mooring order strategically determined using Q in compared to the mooring order strategically determined using Q out , S in , or S out (Figure 9a). As the number of moorings increased, the magnitude of Q in,out and S in,out converged toward the TEF IJ results. In most cases the rate of convergence was qualitatively similar regardless of which TEF parameter was used to determine the mooring order, with some variation in the point at which convergence within 10% of TEF IJ was reached (Figure 9). One exception that stood out was that the number of moorings (μ) at CR M before Q in converged to within 10% of the TEF IJ value was more than double the minimum number of moorings (μ) to reach the same threshold when the mooring placement was optimized for Q out , S in , and S out (Figure 9e). It is unclear why the mooring placement for Q in was particularly inefficient for capturing the magnitude of Q in , but individual components of the exchange flow and salt flux can exhibit different spatial patterns and timescales of variability (e.g., Lerczak et al., 2006), which likely contributed to a misalignment between the locations with the most variance and those with the greatest magnitude in the exchange flow.

Case Studies: Hypothetical Mooring Deployments, TEF case
Both evenly and strategically distributed mooring approaches tested here are skewed toward numerical modeling applications because of the sampling resolution, which would require many sensors and mooring lines that extend across nearly the full water column. In particular, strategic mooring placement requires extensive a priori knowledge of the system. In an effort to connect this analysis more closely to observations, a specific case study, TEF case , was also examined as described in Section 2.4.4. Note. Convergence is defined here by Q in,out and the freshwater fraction equivalent of S in,out consistently reaching within 10% of the full model value.
These estimates of M are conservative and in some cases M can be lower such as with a greater number of sample depths N (as indicated in each footnote). a For N ≥ 5, S in also converges for M = 1, but not M = 2. b For N ≥ 6, S out also converges for M = 1. c For N ≥ 4, Q in converges for M = 3, but not M = 4. d Q out also converges with as few as M = 4 moorings for N ≥ 7. e For N ≥ 3, Q in also converges for M = 5, but not M = 6. f S in also converges for M = 3. g S out also converges for M = 1, but not M = 2. Figure 7. The correlation coefficient between time series of parameters calculated from the full model resolution TEF IJ compared to strategically placed moorings TEF μJ (black) and evenly distributed moorings TEF MN (gray; N = 10) as the number of moorings is increased.

Table 3 For Each Parameter (M, N) Are the Minimum Number of Evenly Distributed Moorings (M) and the Minimum Number of Depth Samples (N) for Which the Magnitude of TEF MN Parameters Converge Toward the TEF IJ Values Calculated at the Full Model Resolution
The salinity interpolation impacted both the calculated exchange volume transport and salinities ( Figure 10). Similar to the other sampling methods presented, as the number of moorings M increased, TEF case converged to a similar value for M ≥ 4, although not necessarily toward TEF IJ . At SDB M , where salinity stratification is weak, TEF case converged within ≈10% of TEF IJ with ≥4 moorings. In other words, the vertical salinity interpolation method did not matter. At sections SJF EH1 and CR M , however, the results were mixed with the two-layer approximation generally performing better than the linear interpolation with varying sensitivity.
At SJF EH1 , linearly interpolating salinity led to underestimating the magnitude of the exchange volume transport (Figures 10a and 10b). This is possible when opposing currents are classified into the same salinity range, which reduces the net exchange magnitude within that salinity class. In this case, a two-layer S approximation led to values of TEF case closer to TEF IJ than a linear interpolation. The results were fairly insensitive to the interface depth between the two uniform salinity layers, except for S out . At CR M , the two-layer S approximation performed slightly better than the linear interpolation (Figures 10c-10h). Whether the two-layer approximation resulted in an improvement of Q in,out over the linear approximation was sensitive to the interface depth. In neither case did the TEF case exchange salinity S in converge to within 10% of the full TEF IJ value.

Discussion
The results of this study suggest that TEF via in situ moorings can be well approximated in many situations with ≤4 moorings across a channel. However, there are some limitations to this analysis that may constrain the generalization of these results. The magnitude and structure of the exchange flow for each estuary may depend on the specific thresholds, time periods, and cross-sections. The sampling approaches explored here were limited to a small number of relatively simple designs and were not adapted for the particular bathymetry or oceanic conditions at each cross-section. Also the examples chosen, while spanning significant parameter space in terms of estuary type and geometry, were not exhaustive. The practical applicability of the results presented in this study are discussed in the following section in the context of these limitations.
The magnitude and structure of the exchange flow for each model may depend on specific thresholds, time periods, and cross-sections chosen for the calculation. The TEF approach requires the number of salinity bins N bins and a threshold for identifying the cut-off between dividing salinities Q thresh to be defined. While the dividing salinity approach (Lorenz et al., 2019) reduces the sensitivity to these choices (here, N bins = 500 and Q thresh = 0.01), in the weakly stratified San Diego Bay which experiences large variability on tidal time scales (e.g., Figure 4f) the number and location of dividing salinities varied depending on threshold choices (e.g., Figure 5i). However the results presented here were summed over the incoming and outgoing layers and were insensitive to these choices, even at section SDB C . At the other cross-sections, the results were insensitive to variations in N bins and Q thresh . Another consideration is that the exchange flow can vary over time, with sampling frequency, and with model grid resolution (Figures 5a-5c) such that the TEF magnitude may be specific to the particular periods examined (in particular annual vs. seasonal time periods). The conclusions drawn from this analysis focus on comparing the relative change in the TEF across methods with each method applied over the same time period and model Note. Convergence is defined here by Q in,out and the freshwater fraction equivalent of S in,out reaching within 10% of the full model value.  Figure 9. TEF μJ parameters calculated for μ = 1-15 moorings strategically placed across the channel, using the full vertical grid resolution at each mooring. Black lines and gray shading indicate the magnitude of each TEF IJ parameter and the ±10% range (calculated using the freshwater fraction for salinities), respectively. Colored lines indicate the parameter used to determine the order of mooring placement by maximizing the correlation between TEF μJ and TEF IJ . Colored dots mark the first point where each TEF μJ parameter converges to within 10% of the corresponding TEF IJ parameter (e.g., Table 4).
cross-section to reduce the dependence on the specific time period and resolution. However, if only part of the year was examined when surface intrusions were absent at SJF EH1 , for example, the optimal lateral spacing may have been impacted.
Variances in the TEF parameter time series were strongly correlated between methods (R > 0.85) even with a single mooring, except at SJF EH1 (Figure 7). The high correlation is likely due to the strong spring-neap and seasonal variability at CR M and SDB M ( Figure 5) that are spatially coherent as opposed to the seasonal spatial  variability caused by surface intrusions at SJF EH1 . That a single mooring location captures most of the TEF variance also suggests that optimizing the correlation may not be a useful method for strategic mooring placement to measure exchange flows. Even at SJF EH1 where the variance in the exchange flow was improved using the TEF μJ method relative to TEF MN (Figures 7a and 7b), the minimum number of moorings for which the TEF μJ magnitude converged to within 10% of TEF IJ was greater than the minimum number of geometrically distributed moorings (Table 3 compared to 4). Also, at CR M , 10 mooring samples were required for Q in using the TEF μJ method to converge to within 10% of TEF IJ result (Table 4). These results suggest that the locations across each section with the most temporal variance are not the same locations where the majority of the transport occurs. That in most cases 4 evenly distributed moorings captured >90% of the exchange flow (Table 3) and most of the flow variance (Figure 7) also suggests that a straightforward sampling plan-as long as there are sufficient number of moorings and depths sampled-is likely to capture both features (variance and magnitude) of the exchange. This is also supported by the observation that there was generally little difference in the results as mooring placement varied (Figure 9).
While the selected cross-sections span a bay, salt-wedge, and fjord estuary type with differing scales and geometric complexity (Geyer & MacCready, 2014), this subset does not comprehensively cover the full range of estuarine shapes, sizes, and dynamics. Nevertheless, given the range investigated here, the similarity of the results was striking. First, using only a single mooring led to an over-estimate of Q in in every case. This over-estimate could be several times the magnitude of Q in from the full model resolution. In particular for the transects sampled in SJF and SDB, with relatively simple U-shaped channels, using only two laterally distributed moorings that did not sample the channel center led to an underestimate of Q in (e.g., Figures 6a and 6i). While the greatest inflow tended to be concentrated toward the central and deeper channel, the lateral distribution of the outflow was more variable (Figure 3). Second, it was encouraging-from the perspective of capturing exchange flows in a range of systems with limited a priori knowledge-that using evenly distributed moorings performed as well or better than the strategic sampling strategy (Table 3 compared to Table 4) and also that the results converged toward TEF IJ even as the specific sampling locations were varied (Figure 9). The mean number of evenly distributed lateral mooring locations across each channel to resolve (Q in,out ) to within 10% was M = 4.0 ± 1.8 with N = 3.3 ± 1.1 sample depths evenly distributed across the water column and M = 2.0 ± 1.5, N = 2.3 ± 2.1 to resolve S in,out (Table 3).
Further studies of systems dominated by either horizontal or vertical shear would be needed to assess if the sampling resolution might be related to K e and E k . While K e and E k are comparable across systems, the predictive utility of such parameters may be complicated by channel curvature and variable bathymetry and requires a priori knowledge of the flow. However, the degree of stratification did appear to be important to understand how the vertical and horizontal sampling resolution impacted TEF. At SJF EH1 and CR M , TEF MN converged to within 10% of TEF IJ with ≤4 evenly distributed moorings ( Figure 6) while TEF case did not always converge to the same values as TEF IJ (Figure 10). This suggests that when stratification is important more vertical resolution of S reduces the number of moorings needed across the channel. Also, at these stratified sections the result depends on the vertical interpolation of S. In contrast, at SDB M the number of moorings needed to accurately calculate TEF MN (Figure 6) was similar to TEF case ( Figure 10) and was similar between cases A and B. At SDB M with u and S samples evenly distributed vertically the exchange volume flux only converged for N ≥ 4, while in the case studies the TEF case exchange volume fluxes converged to TEF IJ when S was sampled at only two depths, near-surface and near-bottom.
Given that only three estuaries (6 cross-sections) were examined here and that this study utilized numerical model output, the question remains: how realistic would it be to apply the results of these experiments to observations and in other systems? One general limitation of observational studies is the cost of moorings (anchors, line, floats, etc.) as well as of the individual sensors. The results of this study may be useful to constrain an estimate of the sign and possibly the magnitude of the error for sampling studies with fewer moorings. For example, one could extrapolate that a single mooring centered in a deep part of the channel is likely to overestimate the magnitude of the exchange volume transport. The outcome of the case studies also suggest that the salinity does not have to be sampled at the same resolution as the currents to estimate the exchange flow with relatively high accuracy, although as stratification increases, identifying the best salinity interpolation remains a challenge (Aristizábal & Chant, 2015). In channels with high ship traffic, near-surface measurements can be particularly challenging. However, it may be encouraging that the results here demonstrated relatively little sensitivity to the specific sample location, that is, one could place a mooring outside of a navigational channel.

Conclusions
Exchange between estuaries and the coastal ocean or through inland seas is an important driver of the circulation, mixing, biology, and chemistry on both ends of the exchange. Significant progress has been made in calculating this exchange in estuarine conditions with time-varying stratification and flow, strong vertical and lateral velocity shear and salinity gradients, and complex bathymetry using the TEF method (Lorenz et al., 2019(Lorenz et al., , 2020. However, application of this theory has predominantly utilized numerical models where the salinity and velocity are highly resolved in space and time. In this analysis TEF was calculated using various methods to sub-sample realistic numerical models in order to understand the sensitivity of TEF and to develop recommendations for minimal sampling thresholds that accurately reproduce the exchange flow. Three different estuaries were examined, including San Diego Bay, the Columbia River, and the Salish Sea exchange through the Strait of Juan de Fuca. These examples span a range of estuary types (bay, salt-wedge, and fjord, respectively), scales, depths, and channel bathymetries. Evenly distributed sample locations across the channel, representative of moorings, was the most efficient way to capture the TEF. Three to four moorings were typically the minimum lateral sample distribution required to capture ≥90% of the exchange transport rate Q in and Q out . In most cases, the exchange volume transport was the limiting parameter, requiring more moorings to measure than the exchange flow salinities S in and S out . The minimum vertical resolution to capture ≥90% of the TEF was similar, N ≥ 4, and was also limited by Q in and Q out . Although the exchange calculated by these methods is also dependent on resolving the salinity, the TEF was less sensitive to resolving salinity at the same vertical resolution as velocity and less sensitive to the salinity interpolation method in systems where there was less vertical stratification (e.g., the San Diego Bay, as compared to the Columbia River or the Strait of Juan de Fuca). The TEF could be reproduced by resolving the currents throughout the water column and only sampling salinity near the surface and bottom. In comparison to geometrically distributing moorings across the channel, strategic sampling based on capturing the temporal exchange flow variance did not improve the ability to capture the exchange flow magnitude, likely a result of the fact that locations of strongest exchange flow are often not the locations with the highest variance. This method also requires a priori knowledge of the flow field, and is ambiguous depending on which aspect of the TEF (i.e., Q in,out or S in,out ) is used to calculate the variance and is therefore not a recommended approach for estuary sampling methodology. Overall the results presented here are promising suggesting that TEF can be captured well with a reasonable number of cross-sectionally distributed moorings and sampling depths and can be used to estimate the sign and magnitude of errors. Future work to examine the exchange flow through a wider range and greater number of estuaries and channels could be useful to further refine the number of moorings and the approach to determine the best salinity interpolation.

Data Availability Statement
Model cross-sections and analysis code used in this analysis are hosted at the UCSD Library Digital Collections (https://doi.org/10.6075/J0DJ5FS8).